{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><h1> Markov models</h1></center>\n",
    "\n",
    "# 1. Basic\n",
    "Suppose we have sequences of observations,$X_1,X_2,\\ldots,X_T$. The Markov models assume that $X_t$ captures all the relevant information for predicting the future. Then we can write the joint distribution as follows\n",
    "\\begin{align}\n",
    "p(X_1,X_2,\\ldots,X_T) &=p(X_1)p(X_2|X_1) \\ldots P(X_T|X_{T-1})  \\\\\n",
    "                      &=p(X_1)\\prod_{t=2}^T p(X_t|X_{t-1})\n",
    "\\end{align}                 \n",
    "If we assume the transition function $p(X_t|X_{t-1})$ is independent of time, then the model is called **homogeneous, stationary or time-invariant**. If we assume that the observed variables are discrete, so $X_t \\in \\lbrace 1,\\ldots, K \\rbrace$, this is called a discrete-state of finite-state Markov model.\n",
    "\n",
    "When the $X_t \\in \\lbrace 1,\\ldots, K \\rbrace $, the conditional distribution $p(X_{t+1}|X_t)$ can be written as a $K \\times K $ matrix, which known as the **transition matrix** $\\mathbf{A}$, where $A_{ij}=p(X_{t+1}=j|X_t=i)$, which satisfies $\\mathbf{A}_{ij} \\geq 0,\\, \\sum_j \\mathbf{A}_{ij}=1$. We can also define the $n$ step transition matrix\n",
    "$$\n",
    "\\mathbf{A}_{ij}^{(n)}=p(X_{t+n}=j|X_t=i)\n",
    "$$\n",
    "we can easy to get that\n",
    "$$\n",
    "\\mathbf{A}_{ij}^{(n)}=\\sum_{k} \\mathbf{A}_{ik}^{(n-1)} \\mathbf{A}_{kj}\n",
    "$$\n",
    "which means $\\mathbf{A}^{(n)}=\\mathbf{A}^{(n-1)}\\,\\mathbf{A}$,hence $\\mathbf{A}^{(n)}=\\mathbf{A}^n$\n",
    "\n",
    "Let $\\vec{u}$ be the probability vector which represents the starting distribution. Then the probability that the chain is in after $n$ steps is as follows\n",
    "$$\n",
    "\\vec{u}^{(n)}=\\left(\\mathbf{A}^T\\right)^n \\vec{u}\n",
    "$$\n",
    "\n",
    "# 2. Language modeling\n",
    "The task of language modeling is to estimate the joint probability distributions over sequances of words $X_1,X_2,\\ldots,X_T$. The marginal probabilities $p(X_t=k)=\\pi_k$ are called **unigram statistics**. If we use a first-order Markov model, then $p(X_t=k|X_{t-1}=j)$ is called a **bigram model**. If we use a second-order Markov model, then $p(X_t=k|X_{t-1}=j,X_{t-2}=i)$ is called a  **trigram model**. In general, these are called **n-gram models**.\n",
    "\n",
    "Language models can be used for several things,such as\n",
    "* **Sentence completion**. A language model can predict the next word given the previous words in a sentence.\n",
    "* **Text classification** Any density model can be used as a class-conditional density and hence truned into a (generative) classifier. 0-gram class-conditional density would be equivalent to a naive Bayes classifier.\n",
    "\n",
    "## 2.1 MLE for Markov language models\n",
    "The likelihood function for dataset $D=\\lbrace (x_1,x_2,\\ldots,x_{T_i}) \\rbrace_{i=1}^N $ is as follows\n",
    "\\begin{align}\n",
    "p(D|\\vec{\\pi},\\mathbf{A}) &=\\prod_{i=1}^N \\pi(x_{i1}) \\prod_{t=2}^{T_i} A(x_{i\\,t-1},x_{it}) \\\\\n",
    "                          &=\\prod_{i=1}^N \\left\\{ \\prod_{k=1}^K \\pi_k^{\\mathbb{1}(x_{i1}=k)} \\prod_{t=2}^{T_i} \\prod_{u=1}^K \\prod_{v=1}^K A_{uv}^{\\mathbb{1}(x_{i\\,t-1}=u,x_{it}=v)} \\right\\}\n",
    "\\end{align}\n",
    "Hence the log-likelihood function is\n",
    "$$\n",
    "log\\,p(D|\\vec{\\pi},\\mathbf{A})=\\sum_{k=1}^K N_k log\\,\\pi_k +\\sum_{u=1}^K \\sum_{v=1}^K N_{uv} log\\,A_{uv}\n",
    "$$\n",
    "where $N_k=\\sum_{i=1}^N \\mathbb{1}(x_{i1}=k)$ and $N_{uv}=\\sum_{i=1}^N \\sum_{t=2}^{T_i} \\mathbb{I}(x_{i\\,t-1}=u,x_{it}=v)$\n",
    "\n",
    "Hence we can write the MLE as the normalizad counts\n",
    "\\begin{align}\n",
    "\\hat{\\pi}_k=\\frac{N_k}{\\sum_{k'=1}^K N_{k'}} \\\\\n",
    "\\hat{\\mathbf{A}}_{jk}=\\frac{N_{jk}}{\\sum_{k'=1}^K N_{jk'}}\n",
    "\\end{align}\n",
    "These first-order result can be extended in a straightforward way to higher order Markov models. However, the problem of zero-counts becomes very acute as the size of vocabulary or the order of the chain becomes large. We need more bayesian.\n",
    "\n",
    "## 2.2 Bayesian version of parameter estimation\n",
    "First, we can use an Dirichlet prior for $\\vec{\\pi}$ and use an independent Dirichlet prior on each row of transition matrix\n",
    "$$\n",
    "\\mathbf{A}_j \\sim Dir(\\alpha_0 \\vec{m})=Dir(\\vec{\\alpha})\n",
    "$$\n",
    "where $\\vec{m}$ is the prior mean (satisfying $\\sum_k m_k=1$) and $\\alpha_0$ is the prior strength. The posterior is given by $\\mathbf{A}_j \\sim Dir(\\vec{\\alpha}+\\vec{N}_j)$,where $\\vec{N}_j=(N_{j1},\\ldots,N_{jK})$ is the vector which records the number of times we have transitioned out of state $j$ to each of the other states. The posterior predictive density is\n",
    "$$\n",
    "p(X_{t+1}=k|X_t=j,D)=\\frac{N_{jk}+\\alpha_0 m_k}{N_j+\\alpha_0}\n",
    "$$\n",
    "The only remaining question is: what value should we use for $\\vec{\\alpha}$. Let's use the empirical Bayes.\n",
    "\\begin{align}\n",
    "p(D|\\vec{\\alpha})=\\int p(D|\\vec{\\pi},\\mathbf{A})p(\\mathbf{A}|\\vec{\\alpha})d\\vec{\\alpha}\n",
    "\\end{align}\n",
    "we can maximum $p(D|\\vec{\\alpha})$ to get the optimal $\\vec{\\alpha}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The true transition matrix is:\n",
      "[[0.12619616 0.32033016 0.24663533 0.20170978 0.05256834 0.05256022]\n",
      " [0.01801633 0.26866992 0.18645344 0.21962946 0.00638488 0.30084597]\n",
      " [0.37179051 0.09483616 0.08120775 0.08191322 0.13588249 0.23436986]\n",
      " [0.20250327 0.13653324 0.28684718 0.06539713 0.13696244 0.17175674]\n",
      " [0.17581597 0.30268703 0.07697468 0.19823849 0.2283771  0.01790673]\n",
      " [0.17036983 0.04781896 0.01824199 0.26608976 0.27078587 0.22669358]]\n",
      "The generated Markov chian is :\n",
      "[6 6 6 ... 4 5 4]\n",
      "The estimated transition matrix is (MLE):\n",
      "[[0.10559006 0.33540373 0.24223602 0.17391304 0.05590062 0.08695652]\n",
      " [0.02439024 0.27317073 0.16585366 0.24878049 0.00487805 0.28292683]\n",
      " [0.37671233 0.12328767 0.08219178 0.08219178 0.13013699 0.20547945]\n",
      " [0.18343195 0.13017751 0.28994083 0.0591716  0.13017751 0.20710059]\n",
      " [0.18518519 0.31111111 0.08148148 0.19259259 0.20740741 0.02222222]\n",
      " [0.15300546 0.07103825 0.00546448 0.2295082  0.31147541 0.2295082 ]]\n"
     ]
    }
   ],
   "source": [
    "# We will make toy markov chain here, which has 6 states corresponding to the six surface of a dice\n",
    "import numpy as np\n",
    "from pprint import pprint\n",
    "\n",
    "np.random.seed(42)\n",
    "transition_matrix=np.random.rand(6,6)\n",
    "transition_matrix=transition_matrix/np.sum(transition_matrix,axis=1,keepdims=True)\n",
    "print(\"The true transition matrix is:\")\n",
    "print(transition_matrix)\n",
    "\n",
    "def random_categorical(p):\n",
    "    return np.array(np.random.multinomial(1,p)).argmax()+1\n",
    "\n",
    "states=np.array(random_categorical(np.ones(6)*1.0/6)).reshape(-1)  # random choose a init state\n",
    "chain_length=1000\n",
    "for i in range(chain_length):\n",
    "    states=np.append(states,random_categorical(transition_matrix[states[-1]-1]))\n",
    "print(\"The generated Markov chian is :\")\n",
    "print(states)\n",
    "\n",
    "#Estimate the transition matrix(MLE)\n",
    "estimated_transition_matrix=np.zeros((6,6))\n",
    "for i in range(states.size-2):\n",
    "    estimated_transition_matrix[states[i]-1,states[i+1]-1]+=1\n",
    "estimated_transition_matrix=estimated_transition_matrix/ np.sum(estimated_transition_matrix,\n",
    "                                                                axis=1,keepdims=True)\n",
    "print(\"The estimated transition matrix is (MLE):\")\n",
    "print(estimated_transition_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Hidden Markov models\n",
    "A **hidden Markov model** or **HMM** consists of a discrete-time, discrete-state Markov chain, with hidden states $z_t \\in \\lbrace 1,\\ldots,K \\rbrace$, plus an observation model $p(x_t|z_t)$. The corresponding joint distribution has the form\n",
    "\\begin{align}\n",
    "p(z_{1:T},x_{1:T}) &= p(z_1)p(x_1|z_1) p(z_2|z_1)p(x_2|z_2) \\ldots p(z_T|z_{T-1})p(x_T|z_T) \\\\\n",
    "                   &=\\left[ p(z_1) \\prod_{t=2}^T p(z_t|z_{t-1}) \\right] \\left[\\prod_{t=1}^T p(x_t|z_t) \\right]\n",
    "\\end{align}\n",
    "The observations $x_t$ in an HMM can be discrete or continuous.\n",
    "\n",
    "## 3.1 Inference in HMMs\n",
    "We now discuss how to infer the hidden state sequence of an HMM, assuming the parameters are known. There are several different kinds of inferential tasks for an HMM.\n",
    "* **Filtering** means to compute the **belief state** $p(z_t|x_{1:t})$ online, as the data streams in. This is called \"filtering\" because it reduces the noise more than simply estimating the hidden state using the current estimate, $p(z_t|x_t)$.\n",
    "* **Smoothing** means to compute $p(z_t|x_{1:T})$ offline. By conditioning on past and future data, our uncertainty will be significantly reduced.\n",
    "* **Fixed lag smoothing** is an interesting compromise between online and offline estimation; it involves computing $p(z_{t-l}|x_{1:t})$ , where $l>0$ called the lag. This gives better performance than filtering, but incurs a slight delay.\n",
    "* **Prediction** Instead of predicting the past given the future, as in fixed lag smoothing, we might want to predict the future given the past, i.e., to compute $p(z_{t+h}|x_{1:t})$,where $h>0$ is called the prediction horizon. For example, suppose $h=2$ we have\n",
    "$$\n",
    "p(z_{t+2}|x_{1:t})=\\sum_{z_{t+1}}\\sum_{z_t} p(z_{t+2}|z_{t+1}) p(z_{t+1}|z_t)p(z_t|x_{1:t})\n",
    "$$\n",
    "The quantity $p(z_{t+h}|x_{1:t})$ is a prediction about future hidden states; it can be converted into a prediction about future observations using\n",
    "$$\n",
    "p(x_{t+h}|x_{1:t})=\\sum_{z_{t+h}} p(x_{t+h}|z_{t+h})p(z_{t+h}|x_{1:t})\n",
    "$$\n",
    "* **MAP estimation** This means computing $\\underset{z_{1:T}}{argmax}\\,p(z_{1:T}|x_{1:T})$ which is a most probable state sequence.\n",
    "* **Posterior samples** If there is more than one plausible interpretation of the data, it can be useful to sample from the posterior $p(z_{1:T}|x_{1:T})$, which contain much more information than the sequence of marginals computed by smoothing.\n",
    "* **Probability of the evidence** We can compute the **probability of the evidence**, $p(x_{1:T})$ by summing up over all hidden paths, $p(x_{1:T})=\\sum_{z_{1:T}} p(z_{1:T},x_{1:T})$, which can be used to classify sequences.\n",
    "\n",
    "### 3.2 The forwards algorithm\n",
    "We now describe how to **recursively** compute the filtered marginals, $p(z_t|x_{1:t})$ in an HMM, which means\n",
    "\\begin{align}\n",
    "p(z_t|x_{1:t})  &=p(z_t|x_{1:t-1},x_t)  \\\\\n",
    "                &=\\frac{p(z_t,x_t|x_{1:t-1})}{p(x_t)} \\\\\n",
    "                &=\\frac{p(z_t|x_{1:t-1})p(x_t|z_t,x_{1:t-1})}{p(x_t)}  \\\\\n",
    "                &=\\frac{p(z_t|x_{1:t-1})p(x_t|z_t)}{p(x_t)}\n",
    "\\end{align}\n",
    "The $p(z_t|x_{1:t-1})$ can be viewed as a prior distribution for $z_t$.The algorithm has two steps. First comes the prediction step, in which we compute the **one-step-ahead predictive density**.\n",
    "$$\n",
    "p(z_t=j|x_{1:t-1})=\\sum_{z_{t-1}} p(z_t=j|z_{t-1})p(z_{t-1}|x_{1:t-1})\n",
    "$$\n",
    "Next comes the update step,in which we compute the $p(z_t|x_{1:t})$ using\n",
    "\\begin{align}\n",
    "p(z_t|x_{1:t}) &=\\frac{p(z_t|x_{1:t-1})p(x_t|z_t)}{p(x_t)} \\\\\n",
    "               &=\\frac{p(z_t|x_{1:t-1})p(x_t|z_t)}{\\sum_{z_t} p(x_t|z_t)p(z_t|x_{1:t-1})}\n",
    "\\end{align}\n",
    "The distribution $p(z_t|x_{1:t})$ is called the (filtered) **belief state** at time t.\n",
    "\n",
    "In addition to computing the hidden states, we can use this algorithm to compute the log probability of the evidence\n",
    "\\begin{align}\n",
    "log\\,p(x_{1:T}|\\theta) &=log\\,p(x_1|\\theta)+log\\, p(x_2|x_1,\\theta) \\ldots+log\\, p(x_T|x_{1:T-1},\\theta) \\\\\n",
    "                       &=log\\,\\left(\\sum_{z_1} p(x_1|z_1,\\theta)p(z_1|\\theta) \\right) +\\sum_{t=2}^T log\\,\\left(\\sum_{z_t}p(x_t|z_t)p(z_t|x_{1:t-1}) \\right)\n",
    "\\end{align}\n",
    "\n",
    "## 3.3 The forwards-backwards algorithm\n",
    "We now discuss how to compute the smoothed marginals,$p(z_t=j|x_{1:T})$ using offline inference. The basic idea is as follows\n",
    "\\begin{align}\n",
    "p(z_t|x_{1:T}) &=p(z_t|x_{1:t},x_{t+1:T}) \\\\\n",
    "               &=\\frac{p(z_t,x_{t+1:T}|x_{1:t})}{p(x_{t+1:T}|x_{1:t})} \\\\\n",
    "               &=\\frac{p(z_t|x_{1:t})p(x_{t+1:T}|z_t,x_{1:t})}{p(x_{t+1:T}|x_{1:t})} \\\\\n",
    "               &=\\frac{p(z_t|x_{1:t})p(x_{t+1:T}|z_t)}{p(x_{t+1:T}|x_{1:t})} \\\\\n",
    "\\end{align}\n",
    "We have using the forwards algorithm to compute the $p(z_t|x_{1:t})$.The conditional likelihood of future evidence given that the hidden state at time t is $p(x_{t+1:T}|z_t)$,which can also compute using **recursive**, but in a backwards manner.\n",
    "\\begin{align}\n",
    "p(x_{t+1:T}|z_t) &=\\sum_{z_{t+1}} p(x_{t+1:T}|z_{t+1})p(z_{t+1}|z_t) \\\\\n",
    "                 &=\\sum_{z_{t+1}} p(x_{t+2:T},x_{t+1}|z_{t+1})p(z_{t+1}|z_t) \\\\\n",
    "                 &=\\sum_{z_{t+1}} p(x_{t+2:T}|z_{t+1})p(x_{t+1}|z_{t+1}) p(z_{t+1}|z_t) \\\\\n",
    "\\end{align}\n",
    "\n",
    "### Two-slice smoothed marginals\n",
    "When we estimate the parameters of the transition matrix using EM,we will need to compute the expected number of transitions from state $i$ to state$j$.\n",
    "\\begin{align}\n",
    "N_{ij} &=\\sum_{t=1}^{T-1} \\mathbb{E} \\left[\\mathbb{I}(z_t=i,z_{t+1}=j)|x_{1:T} \\right] \\\\\n",
    "       &=\\sum_{t=1}^{T-1} p(z_t=i,z_{t+1}=j|x_{1:T})\n",
    "\\end{align}\n",
    "The term $p(z_t=i,z_{t+1}=j|x_{1:T})$ is called a smoothed **two-slice marginal**, and can be computed as \n",
    "\\begin{align}\n",
    "p(z_t=i,z_{t+1}=j|x_{1:T})&=p(z_t=i|x_{1:T})p(z_{t+1}=j|z_t=i,x_{1:T}) \\\\\n",
    "                          &=p(z_t=i|x_{1:T})p(z_{t+1}=j|z_t=i,x_{t+1:T}) \\\\\n",
    "                          &=p(z_t=i|x_{1:T})\\frac{p(z_{t+1}=j|z_t=i)p(x_{t+1:T}|z_{t+1}=j,z_t=i)}{p(x_{t+1:T}|z_t=i)}\\\\\n",
    "                          &=p(z_t=i|x_{1:T})\\frac{p(z_{t+1}=j|z_t=i)p(x_{t+1:T}|z_{t+1}=j)}{p(x_{t+1:T}|z_t=i)}\\\\\n",
    "                          &=p(z_t=i|x_{1:T})\\frac{p(z_{t+1}=j|z_t=i)p(x_{t+1:T}|z_{t+1}=j)}{\\sum_{z_{t+1}} p(z_{t+1}=j|z_t=i)p(x_{t+1:T}|z_{t+1}=j)}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The true transition matrix between three dices is:\n",
      "[[0.18205878 0.46212909 0.35581214]\n",
      " [0.65738127 0.17132261 0.17129612]\n",
      " [0.03807826 0.56784481 0.39407693]]\n",
      "The true observation matrix is:\n",
      "[[0.24206172 0.00703702 0.33157343 0.28457888 0.07259026 0.06215869]\n",
      " [0.07812991 0.12960651 0.22354506 0.18400761 0.12406296 0.26064796]\n",
      " [0.06230408 0.13048463 0.16363328 0.20370088 0.35069405 0.08918308]]\n",
      "The forward algorithm estimate hidden state(MAP) accuracy : 0.537800\n",
      "The forward-backwards algorithm estimate hidden state(MAP) accuracy : 0.556200\n"
     ]
    }
   ],
   "source": [
    "# We will make toy hidden markov chain here, which has 3 hidden states corresponding to the 3 different dice\n",
    "import numpy as np\n",
    "from pprint import pprint\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "n_hidden_states=3\n",
    "np.random.seed(42)\n",
    "transition_matrix=np.random.rand(n_hidden_states,n_hidden_states)\n",
    "transition_matrix=transition_matrix/np.sum(transition_matrix,axis=1,keepdims=True)\n",
    "print(\"The true transition matrix between three dices is:\")\n",
    "print(transition_matrix)\n",
    "\n",
    "observed_prob=np.random.rand(n_hidden_states,6)\n",
    "observed_prob=observed_prob/ np.sum(observed_prob,axis=1,keepdims=True)\n",
    "print(\"The true observation matrix is:\")\n",
    "print(observed_prob)\n",
    "\n",
    "def random_categorical(p):\n",
    "    return np.array(np.random.multinomial(1,p)).argmax()+1\n",
    "\n",
    "#random choose a init hidden state\n",
    "hidden_states=np.array(random_categorical(np.ones(n_hidden_states)*1.0/n_hidden_states)).reshape(-1)  \n",
    "observaed_variable=np.array(random_categorical(observed_prob[hidden_states[-1]-1])).reshape(-1)\n",
    "\n",
    "chain_length=5000\n",
    "for i in range(chain_length-1):\n",
    "    hidden_states=np.append(hidden_states,random_categorical(transition_matrix[hidden_states[-1]-1]))\n",
    "    observaed_variable=np.append(observaed_variable,random_categorical(observed_prob[hidden_states[-1]-1]))\n",
    "'''print(\"The generated observed variables are:\")\n",
    "print(observaed_variable)\n",
    "print(\"The underline hidden states are:\")\n",
    "print(hidden_states)'''\n",
    "\n",
    "def vec_norm(vec): # l1-norm of a vector\n",
    "    return vec/np.sum(vec,keepdims=True)\n",
    "\n",
    "def calculate_z_i(x_i,observed_prob):  # calculate p(z_i|x_i)\n",
    "    return vec_norm(observed_prob[:,x_i-1])\n",
    "\n",
    "# Forward algorithm\n",
    "hidden_forward_estimate=np.zeros((chain_length,n_hidden_states))\n",
    "hidden_forward_estimate[0,:]=calculate_z_i(observaed_variable[0],observed_prob)\n",
    "for step in np.arange(1,chain_length):\n",
    "    un_normalized=[]\n",
    "    for h in np.arange(1,n_hidden_states+1):\n",
    "        un_normalized.append(observed_prob[h-1,observaed_variable[step]-1]*\n",
    "                             np.dot(hidden_forward_estimate[step-1,:],transition_matrix[:,h-1]))\n",
    "    hidden_forward_estimate[step,:]=vec_norm(np.array(un_normalized))\n",
    "    \n",
    "acc=accuracy_score(hidden_forward_estimate.argmax(axis=1)+1,hidden_states)\n",
    "print(\"The forward algorithm estimate hidden state(MAP) accuracy : %.6f\" %acc)\n",
    "\n",
    "# Forward-backwards algorithm\n",
    "backwards_evidence=np.ones((chain_length,n_hidden_states))\n",
    "for step in np.arange(0,chain_length-1)[::-1]:\n",
    "    un_normalized=[]\n",
    "    for h in np.arange(1,n_hidden_states+1):\n",
    "        un_normalized.append(np.sum(backwards_evidence[step+1]* \n",
    "                                    observed_prob[:,observaed_variable[step+1]-1]* transition_matrix[h-1,:]))\n",
    "    backwards_evidence[step,:]=vec_norm(np.array(un_normalized))\n",
    "forward_backwards_estimate=hidden_forward_estimate*backwards_evidence\n",
    "acc=accuracy_score(forward_backwards_estimate.argmax(axis=1)+1,hidden_states)\n",
    "print(\"The forward-backwards algorithm estimate hidden state(MAP) accuracy : %.6f\" %acc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 The Viterbi algorithm\n",
    "The Viterbi algorithm can be used to compute the most probable sequence of states in a chain-structured graphical model.\n",
    "\\begin{align}\n",
    "\\hat{\\vec{z}} &=\\underset{z_{1:T}}{argmax}\\,p(z_{1:T}|x_{1:T})   \\\\\n",
    "\\end{align}\n",
    "This is a **dynamic programming** problem. We have to find the optimal subproblem.\n",
    "\\begin{align}\n",
    "\\delta_t(j) &=\\underset{z_1,z_2,\\ldots,z_{t-1}}{max}\\,p(z_{1:t-1},z_t=j|x_{1:t})  \\\\\n",
    "            &=\\underset{z_1,z_2,\\ldots,z_{t-1}}{max}\\,p(z_{1:t-1}|x_{1:t})p(z_t=j|z_{t-1},x_t) \\\\\n",
    "            &=\\underset{z_1,z_2,\\ldots,z_{t-1}}{max}\\,p(z_{1:t-1}|x_{1:t})\\frac{p(x_t|z_t=j)p(z_t=j|z_{t-1})}{p(x_t|z_{t-1})}\\\\\n",
    "            &=\\underset{z_1,z_2,\\ldots,z_{t-1}}{max}\\,p(z_{1:t-1}|x_{1:t-1},x_t)\\frac{p(x_t|z_t=j)p(z_t=j|z_{t-1})}{p(x_t|z_{t-1})}\\\\\n",
    "            &=\\underset{z_1,z_2,\\ldots,z_{t-1}}{max}\\,\\frac{p(x_t|z_{t-1})p(z_{1:t-1}|x_{1:t-1})}{p(x_t|x_{1:t-1})} \\frac{p(x_t|z_t=j)p(z_t=j|z_{t-1})}{p(x_t|z_{t-1})}\\\\\n",
    "            &=\\underset{i}{max} \\underset{z_1,z_2,\\ldots,z_{t-2}}{max}\\,\\frac{p(x_t|z_{t-1}=i)p(z_{1:t-2},z_{t-1}=i|x_{1:t-1})}{p(x_t|x_{1:t-1})} \\frac{p(x_t|z_t=j)p(z_t=j|z_{t-1}=i)}{p(x_t|z_{t-1}=i)}\\\\\n",
    "            &=\\underset{i}{max} \\underset{z_1,z_2,\\ldots,z_{t-2}}{max}\\,\\frac{p(x_t|z_{t-1}=i)\\delta_{t-1}(i)}{p(x_t|x_{1:t-1})} \\frac{p(x_t|z_t=j)p(z_t=j|z_{t-1}=i)}{p(x_t|z_{t-1}=i)}\\\\\n",
    "            &=\\underset{i}{max} \\underset{z_1,z_2,\\ldots,z_{t-2}}{max}\\,\\frac{\\delta_{t-1}(i)p(x_t|z_t=j)p(z_t=j|z_{t-1}=i)}{p(x_t|x_{1:t-1})}\\\\\n",
    "            &=\\underset{i}{max} \\underset{z_1,z_2,\\ldots,z_{t-2}}{max}\\,\\delta_{t-1}(i)p(x_t|z_t=j)p(z_t=j|z_{t-1}=i)\n",
    "\\end{align}\n",
    "So the optimal subproblem is $\\delta_t(j)=\\underset{z_1,z_2,\\ldots,z_{t-1}}{max}\\,p(z_{1:t-1},z_t=j|x_{1:t})$, we can easily use the dynamic programming technology to solve this problem.\n",
    "\n",
    "### Forwards filtering, backwards sampling\n",
    "It is often useful to sample paths from the posterior\n",
    "$$\n",
    "z_{1:T}^s \\sim p(z_{1:T}|x_{1:T})\n",
    "$$\n",
    "The key insight into how to do this is that we can write the joint from right to left using\n",
    "$$\n",
    "p(z_{1:T}|x_{1:T})=p(z_T|x_{1:T}) \\prod_{t=T-1}^1 p\\left(z_t|z_{t+1},x_{1:t}\\right)\n",
    "$$\n",
    "We can then sample $z_t$ given future sampled states using\n",
    "$$\n",
    "z_t^s \\sim p(z_t|z_{t+1}^s,x_{1:t})=\\frac{p(z_t|x_{1:t})p(z_{t+1}^s|z_t)}{p(z_{t+1}^s|x_{1:t})}\n",
    "$$\n",
    "we need to do forwards filtering first.\n",
    "\n",
    "# 4. Learning for HMMs\n",
    "We now discuss how to estimate the parameter $\\mathbf{\\theta}=\\left(\\pi,\\mathbf{A},\\mathbf{B} \\right)$, where $\\pi_i=p(z_1=i)$ is the initial state distribution, $A(i,j)=p(z_t=j,z|z_{t-1}=i)$ is the transition matrix, and $\\mathbf{B}$ are the parameters of the conditional densites $p(x_t|z_t=j)$.\n",
    "\n",
    "## 4.1 Training with fully observed data\n",
    "If we observe the hidden state sequences, we can compute the MLEs for $\\mathbf{A}$ and $\\pi$ exactly as in the Markov models. The details on how to estimate $\\mathbf{B}$ depend on form of the observation model which is identical to fitting a generative classifier. \n",
    "\n",
    "## 4.2 EM for HMMs (the Baum-Welch algorithm)\n",
    "If the $z_t$ variable are not observed, we are in a situation analogous to fitting a mixture model.The most common approach is to use the EM algorithm to find the MLE or MAP parameters. \n",
    "#### E step\n",
    "The complete data log likelihood function is as follows\n",
    "\\begin{align}\n",
    "log\\,p(D|\\mathbf{\\theta}) &=\\sum_{i=1}^N \\left[ log\\,p(z_{i1})+\\sum_{t=1}^{T_i} log\\,p(x_{it}|z_{it}) +\\sum_{t=2}^{T_{t}} log\\, p(z_{it}|z_{i,t-1}) \\right] \\\\\n",
    "                          &=\\sum_{i=1}^N \\left[ log\\,\\prod_{k=1}^K \\pi_k^{\\mathbb{1}(z_{i1}=k)} +\\sum_{t=1}^{T_i} log\\,p(x_{it}|z_{it})+\\sum_{t=2}^{T_i} log\\,\\prod_{u=1}^k\\prod_{v=1}^k A_{uv}^{\\mathbb{1}(x_{i,t-1}=u,x_{i,t}=v)} \\right] \\\\\n",
    "\\end{align}\n",
    "we can get the expected complete data log likelihood as\n",
    "$$\n",
    "\\mathcal{Q}(\\theta,\\theta^{old}) =\\mathbb{E}\\left[log\\,p(D|\\mathbf{\\theta}) \\right] \n",
    "$$\n",
    "where the expectation is taken respect to $\\vec{z_i} \\sim p(\\vec{z}_i|\\vec{x}_i,\\theta^{old})$, so we can rewrite it as\n",
    "\\begin{align}\n",
    "\\mathcal{Q}(\\theta,\\theta^{old}) &=\\mathbb{E}\\left[log\\,p(D|\\mathbf{\\theta}) \\right] \\\\\n",
    "                                 &=\\mathbb{E}\\left\\{\\sum_{i=1}^N \\left[ log\\,\\prod_{k=1}^K \\pi_k^{\\mathbb{1}(z_{i1}=k)} +\\sum_{t=1}^{T_i} log\\,p(x_{it}|z_{it})+\\sum_{t=2}^{T_i} log\\,\\prod_{u=1}^k\\prod_{v=1}^k A_{uv}^{\\mathbb{1}(x_{i,t-1}=u,x_{i,t}=v)} \\right] \\right\\} \\\\\n",
    "                                 &=\\sum_{k=1}^K \\mathbb{E}\\left[N_k^1\\right]log\\,\\pi_k+\\sum_{i=1}^N\\sum_{t=1}^{T_i}\\sum_{k=1}^K p(z_{it}=k|\\vec{x}_i,\\theta_{old})log\\,p(x_{it}|\\phi_k)+\\sum_{j=1}^K\\sum_{k=1}^K\\mathbb{E}\\left[N_{jk}\\right]log\\,A_{jk} \\\\\n",
    "\\end{align}\n",
    "where we have defined\n",
    "\\begin{align}\n",
    "\\mathbb{E}\\left[N_k^1\\right] &=\\sum_{i=1}^N p(z_{i1}=k|\\vec{x}_i,\\theta_{old}) \\\\\n",
    "\\mathbb{E}\\left[N_{jk}\\right] &=\\sum_{i=1}^N \\sum_{t=2}^{T_i}p(z_{i,t-1}=j,z_{it}=k|\\vec{x}_i,\\theta_{old})\n",
    "\\end{align}\n",
    "which can be calculated using forwards-backwards algorithm.\n",
    "#### M step\n",
    "It's easy to get the maximum estimate for $\\mathbf{A}$ and $\\pi$ as follows\n",
    "\\begin{align}\n",
    "\\hat{\\mathbf{A}}_{jk} &=\\frac{\\mathbb{E}\\left[N_{jk}\\right]}{\\sum_{k'} \\mathbb{E}\\left[N_{jk'}\\right]} \\\\\n",
    "\\hat{\\pi}_k           &=\\frac{\\mathbb{E}\\left[N_k^1\\right]}{N} \n",
    "\\end{align}\n",
    "The maxmium estimate for $\\phi_k$ is depended on the form of the conditional distribution $p(x_{it}|\\phi_k)$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
